The metaTF package is a computational framework designed to perform transcription factor (TF) activity analysis using single-cell RNA sequencing data. This vignette introduced the basic usage of metaTF for our test data. For more usage of this package, please refer to the following vignette:

Case study I: Identifying cell-type specific TF regulons during HSC formation

Case study II: Decoding the dynamic TF activity during HSC formation

1. Quick start

Here we list the simplified usage of metaTF:

suppressMessages(suppressWarnings(library(metaTF)))
# 1. load data
exp_data <- readRDS(system.file("./extdata/test_exp_data.rds",package = "metaTF"))
col_data <- readRDS(system.file("./extdata/test_col_data.rds",package = "metaTF"))
# 2. create object and inferring GRNs
atfr <- metaTF(exp_data = exp_data, col_data = col_data)
atfr <- scater::logNormCounts(atfr)
atfr <- inferGRNs(x = atfr, method="PUIC", ncores=4)
# 3. integrating GRNs and get regulons
PUIC_net <- PUIC::matToNet(weightMat = as.matrix(GRN(x = atfr,i = "PUIC")), methods = "aracne") %>% filter(regulator %in% TF_Symbols$mouse)
## [1] Generating Gene regulatory networks using: aracne.
blood_net <- readRDS(system.file("./extdata/mm_blood_network.rds",package = "metaTF")) %>% filter(tf %in% PUIC_net$regulator)
merged_grns <- integraNets(grnNets = list(PUIC_net, blood_net), norm_weight = T, return_sig = T, ncores = 6, max_target_num = 200, maxInter = 10)
merged_grns <- dfToList(merged_grns)
atfr <- filterRegulons(x = atfr, gene_list = merged_grns, use_grn="PUIC")
## Filtering regulons for 61 TFs.
# 4. regulon activity estimation
atfr <- regulonActivity(x = atfr)
## 
## Computing the association scores
## Computing regulons enrichment with aREA
# 5. dimension reduction
atfr <- regulonPCA(x = atfr)
## Warning in (function (A, nv = 5, nu = nv, maxit = 1000, work = nv + 7, reorth =
## TRUE, : You're computing too large a percentage of total singular values, use a
## standard svd instead.
atfr <- regulonTSNE(x = atfr)
atfr <- regulonUMAP(x = atfr)
p <- plotRegulonReducedDim(object = atfr, alt_assay = "viper", dimred = "PCA", colour_by="stage")
q <- plotRegulonReducedDim(object = atfr, alt_assay = "viper", dimred = "TSNE", colour_by="stage")
r <- plotRegulonReducedDim(object = atfr, alt_assay = "viper", dimred = "UMAP", colour_by="stage")
p+q+r

# 6. cell-type specific regulon
cst_regulon <- findAllCTSRegulons(x = atfr, clust_name="stage")
head(cst_regulon)
# 7. regulon-pathway similarity
mm_pathway <- Reactome_pathway$`Mus musculus`
regulon_pathway <- calcuModuleSimilar(tfr_list = as.list(Regulon(x = atfr)), pathway_list = mm_pathway, ncores = 2)
head(regulon_pathway)

2. Creating object

The metaTF package were developed base on S4 Class SingleCellExperiment. Supporting various type of functions which could manipulate the SingleCellExperiment object. We supported the input data such as expression matrix or objects such as Seurat.

2.1 Installing packages and creating object

Here we created an object named atfr firstly using part of the test data from our paper. Our single-cell Smart-seq2 sequencing data for this vignette including 1,000 genes and 154 cells from 6 HSC formation stages.

Firstly, we load expression matrix and a metadata from our package:

#install.packages("devtools", repos = "https://cran.rstudio.com/")
#devtools::install_github("wangdongsum/metaTF")
suppressMessages(suppressWarnings(library(metaTF)))
#importing local expression data
exp_data <- readRDS(system.file("./extdata/test_exp_data.rds",package = "metaTF"))
exp_data[1:4,1:4]
##          E10.0_AEC_1_1 E10.0_AEC_1_2 E10.0_AEC_1_3 E10.0_AEC_1_4
## Nanos1               0             0             0             0
## Lamtor4            102             9             0             9
## Arhgap39             0             0             0             0
## Gm44610              0             0             3             0

We suggest the row names of expression data is gene symbol since the TF regulons are readable in following steps. Next, we import the column annotation data, with each row represent a cell in exp_data matrix. The row names of col_data must be column names of exp_data.

col_data <- readRDS(system.file("./extdata/test_col_data.rds",package = "metaTF"))
head(col_data)

Next, we create an object named atfr using the expression and column metadata.

atfr <- metaTF(exp_data = exp_data, col_data = col_data)
atfr
## class: SingleCellExperiment 
## dim: 1000 154 
## metadata(2): metaTF atfr_version
## assays(1): counts
## rownames(1000): Nanos1 Lamtor4 ... Gm42583 Atp6v0a1
## rowData names(0):
## colnames(154): E10.0_AEC_1_1 E10.0_AEC_1_2 ... Adult_HSC_3739
##   Adult_HSC_3740
## colData names(2): sample stage
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

The metaTF package rely on SingleCellExperiment class, the users could also create object using SingleCellExperiment(), for example:

atfr <- SingleCellExperiment(list(counts=exp_data),colData = col_data)
atfr
## class: SingleCellExperiment 
## dim: 1000 154 
## metadata(0):
## assays(1): counts
## rownames(1000): Nanos1 Lamtor4 ... Gm42583 Atp6v0a1
## rowData names(0):
## colnames(154): E10.0_AEC_1_1 E10.0_AEC_1_2 ... Adult_HSC_3739
##   Adult_HSC_3740
## colData names(2): sample stage
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

2.2 Creating object directly from other objects

Most of the users perform single-cell analysis using Seurat package. We also support data import directly from Seurat object. The users could also choose to keep the reduced dimensions or log-scaled counts or not.

seu <- Seurat::CreateSeuratObject(counts = exp_data)
## Warning: The following arguments are not used: row.names
seu@meta.data$orig.ident <- col_data[row.names(seu@meta.data),]$stage
atfr_seu <- importFromSeurat(object = seu, add_reductions = F)

2.3 Normalizing the expression data

We recommend user to preform analysis using normalized expression data.

#BiocManager::install(pkgs = "scater")
atfr <- scater::logNormCounts(atfr)
assayNames(atfr)
## [1] "counts"    "logcounts"

3. Gene regulatory network (GRN) inferring

The metaTF incorporated several gene regulatory network inferring tools, such as PUIC, GENIE3 or PCOR.

3.1 Inferring global gene regulator network

Here we suggesting using our modified version of PUIC, which is based on information entropy theory and we optimization the running time. Meanwhile, we strongly suggest you to use multi-core computers for expression data with over 5,000 features(genes).

atfr <- inferGRNs(x = atfr, use_assay="logcounts", method = "PUIC", ncores = 6)

Generally, we recommend users NOT to set regulators and targets when choosing PUIC method, because the PUIC algorithm calculated weights based on triple gene pairs, absence of regulators or targets will decrease accuracy.

To view all the regulatory network names in your object, using the GRNs function:

GRNames(atfr)
## [1] "PUIC"

To obtain the weighted GRN matrix:

PUIC_mat <- as.matrix(GRN(x = atfr, i = "PUIC"))
PUIC_mat[1:4,1:4]
##          Nanos1 Lamtor4 Arhgap39 Gm44610
## Nanos1    1.000   0.478    0.206   0.949
## Lamtor4   0.478   1.000    0.627   0.324
## Arhgap39  0.206   0.627    1.000   0.367
## Gm44610   0.949   0.324    0.367   1.000

We also support users to calculate the GRN using their favorite methods and then add the weighted GRN matrix into the atfr object using following function:

GRN(x = atfr, i = "PUIC") <- PUIC_mat

Tips: The inferGRNs() method support 4 methods to inferring GRN, when choosing genie3, we strongly suggest you to set the use_regulator and use_target parameter since this methods is time-consuming when you have too much genes.

3.2 Generating GRN-based network

Here we removed the weakest connections between regulators and targets, leaving the strongest connections from GRN using matToNet method from PUIC package.

PUIC_net <- PUIC::matToNet(weightMat = PUIC_mat, methods = "aracne")
## [1] Generating Gene regulatory networks using: aracne.
head(PUIC_net)

The output is a three-column weighted network with weakest edges removed. Since the matToNet() required square matrix as input. what if you have a weighted network which is not a square matrix (eg. TF by target matrix)? Here we suggest you to keep the top ranked targets (eg. top 5%) for each regulator for further analysis.

library(dplyr)
PUIC_net <- reshape2::melt(data = PUIC_mat)
colnames(PUIC_net) <- c("regulator","target","weight")
PUIC_net_top <- PUIC_net %>% dplyr::group_by(regulator) %>% dplyr::top_n(round(ncol(PUIC_mat)*0.05),weight)
head(PUIC_net_top)

4. Network integration

The metaTF performed network integration from multiple sources, such as databases, motif binding and knockout experiments. Here we introduce the network integration step by step.

4.1 importing prior-based GRN with tissue information

The prior-GRN integrate TF-target information from database such as TRRUST, Chip-seq based resources such as jaspar and TF knockout information from KnockTF database. Here we use mouse blood associated network as prior-GRN for our test data. For complete list of 32 tissue-associated TF-target network (available for human and mouse), please download from Zenodo.

mouse_blood <- readRDS(system.file("./extdata/mm_blood_network.rds",package = "metaTF"))
data(TF_Symbols)
mouse_blood_net <- mouse_blood %>% dplyr::filter(tf %in% TF_Symbols$mouse)

The current prior-GRN at Zenodo only have human and mouse data. For other species, you can map the human prior-GRN to other species using getHomologGRN() and getHomologGene() in our package:

human_pantissue <- readRDS(system.file("./extdata/hs_pantissue_network.rds", package = "metaTF"))
rat_pantissue <- getHomologGRN(df = human_pantissue, input_tax = "human", output_tax = "rat")
rat_tf <- getHomologGene(input_gene = TF_Symbols$human,output_tax = "rat")
rat_pantissue_net <- rat_pantissue %>%dplyr::filter(tf %in% rat_tf$Gene_Symbol_rat)

To browser all the supported species, please use:

browserTaxonomy()

4.2 Network integration and regulon generation

In this step, we integrated the prior-GRN (blood network) with GRN inferred from PUIC using integraNets() and generate

#keeping expressed TFs for further analysis
expressed_tfs <- intersect(row.names(atfr),TF_Symbols$mouse)
#keeping TFs as regulator using GRN from PUIC
PUIC_net <- PUIC_net %>% filter(regulator %in% expressed_tfs)
mouse_blood_net <- mouse_blood_net %>% filter(tf %in% expressed_tfs)
#merge GRNs
merged_grns <- integraNets(grnNets = list(PUIC_net, mouse_blood_net), norm_weight = T, return_sig = T, ncores = 6, max_target_num = 200, maxInter = 10)

Next, we transform the network into a list for further regulon activity evaluation.

grn_list <- dfToList(df = merged_grns)

4.3 (optional) Regulon enrichment using weighted GRN data

If you want to find the most enriched targets of each TF regulon, we suggest you to perform enrichment analysis based on the GRN using filterRegulon().

atfr <- filterRegulons(x = atfr, gene_list = grn_list, use_grn= "PUIC",score_type="pos")
## Filtering regulons for 61 TFs.
filter_grns <- as.list(Regulon(atfr))

Note: when you using your own regulons (eg, dorothea) for the following analysis, we strongly suggest you to perform regulon enrichment analysis, this step will find the most enriched targets for each TF regulon based on your single-cell data.

5. TF regulon activity estimation

Here we incorporated three different methods to estimate the TF regulon activities, the viper and ssgsea performed better in our data.

atfr <- regulonActivity(x = atfr, gene_list = filter_grns, method = "viper", use_assay="logcounts", norm_act=T, norm_coe=10)
## 
## Computing the association scores
## Computing regulons enrichment with aREA
active_mat <- assay(altExp(x = atfr, e = "viper"))

6. Further analysis

In this section, we list several functions implemented in metaTF. Including dimension reduction, cell-type specific regulon identification and TF regulon pathway analysis.

6.1 Dimension reduction of TF regulons

In this section, we performed dimension reduction analysis using TF regulon activity score.

6.1.1 PCA analysis:
atfr <- regulonPCA(x = atfr)
## Warning in (function (A, nv = 5, nu = nv, maxit = 1000, work = nv + 7, reorth =
## TRUE, : You're computing too large a percentage of total singular values, use a
## standard svd instead.
plotReducedDim(object = altExp(atfr),dimred = "PCA",colour_by = "stage")

6.1.2 tSNE analysis:
atfr <- regulonTSNE(x = atfr)
plotReducedDim(object = altExp(atfr, e = "viper"),dimred = "TSNE",colour_by = "stage")

##### 6.1.3 UMAP analysis:

atfr <- regulonUMAP(x = atfr)
plotReducedDim(object = altExp(atfr, e = "viper"),dimred = "UMAP",colour_by = "stage")

Tips: The default order of stage in above figure is in alphabetical order, if you want to change the order of orig.ident in your figure, please use reLevel() in our package. The reLevel() function modified the levels of columns in colData(), including the assays in altExps().

new_order <- c("AEC", "HEC", "T1_pre_HSC", "T2_pre_HSC", "FL_HSC", "Adult_HSC")
atfr <- reLevel(x = atfr, colable = "stage", level_order = new_order)
plotReducedDim(altExp(atfr), dimred = "UMAP",colour_by = "stage")

Meanwhile, the reLevel() also support Seurat objects:

seu <- NormalizeData(object = seu)
seu <- FindVariableFeatures(object = seu)
seu <- ScaleData(object = seu)
## Centering and scaling data matrix
seu <- RunPCA(object = seu)
## PC_ 1 
## Positive:  Psmb8, Nckap1l, Slfn4, Pik3cd, Pou2f2, Tnfrsf18, Dennd1c, Srsf5, Fam110a, Coro2a 
##     Tmem40, Socs2, Retreg1, Angpt1, Sema4a, Sorl1, Samd9l, Tmem8, Cd164, Stat6 
##     H2-DMb1, Igtp, Olfr90, Ppt1, Vmn2r114, Anxa11, Mapk13, 9930111J21Rik2, Gprc5b, Gm13613 
## Negative:  Cks1b, Ppia, Plxnd1, Eif4a1, Psmb6, Nme1, Ipo11, Dnmt3a, Eif4g1, Ltbp1 
##     Adgrf5, Eif4a3, Actl6a, Psmb7, Erp29, Piezo1, Abhd17a, Chchd2, Tubb5, Ddx39b 
##     Kdr, Myh10, Dnmt1, Srsf7, Ywhag, Gja4, Ndufa11, Mpzl1, Igf2r, Fkbp10 
## PC_ 2 
## Positive:  Gja4, Oaz1, Rnasek, Kdr, Htra1, Sipa1l2, Sh3bgrl3, Sema3g, Nos3, Pea15a 
##     App, Rps15, Fbln2, Rabac1, Sec61b, Fignl2, Map1b, Sulf1, Gm2000, Cps1 
##     Timm13, Cdc42ep5, Cspg4, Insyn1, Col1a1, Ndufc1, Hras, Serpina1b, Map6, Fam107a 
## Negative:  mt-Atp6, mt-Nd5, Xpo1, Srsf1, Angpt1, Rsbn1l, Nucks1, Cd164, Fam8a1, Ptp4a1 
##     Fbxo11, Atrx, Gm13212, Naa15, Dhx36, Dok2, Vps35, Cstf2t, Pds5b, Cep192 
##     Gm37697, Eif2s3x, Plcg2, Bclaf1, AI506816, Ap1g1, Mis18bp1, Ddx20, Bdp1, Vcpip1 
## PC_ 3 
## Positive:  Nos3, Sipa1l2, Sall4, Sema3g, Plxnd1, Htra1, Kdr, Pea15a, Fignl2, Naxe 
##     Ndufa11, Rnasek, Lama3, Sh3bgrl3, Gja4, Urod, Glg1, Csk, Myh10, Nisch 
##     Ldlrad4, Taok2, Ccdc28b, Scrn2, Laptm4b, Pin1, Bambi, Hras, Abhd17a, Gm2000 
## Negative:  Adcy2, Hcn2, Cyp2a4, Steap4, A1bg, Ttr, Myt1l, Lamp5, Inmt, Rarres2 
##     Ftcd, Mup20, Mup11, Slc27a5, Afm, Ahsg, Cps1, Kng2, Cyp2j5, Serpina1b 
##     Cadm2, Ina, Afp, Onecut2, Col1a1, Prr7, Thrsp, Mgat3, Ggt7, Map6 
## PC_ 4 
## Positive:  Slfn4, Sema4a, Itgal, Trim11, L1cam, Mapk13, Ostf1, Sorl1, Adpgk, Gm44873 
##     Abtb1, Rnf149, Stat6, Iqgap1, Incenp, Mybl2, Mical1, Spag5, Ccna2, Emb 
##     A230085B16Rik, Mis18bp1, Gm49897, Shcbp1, Dnmt1, Arhgap19, Ccnb1, Gm42519, Iqgap3, Cpd 
## Negative:  Rpl13a, Rpl27a, Rps4x, Rps13, Btf3, Hsp90ab1, Sod1, Selenop, Rps15, Ptp4a2 
##     Mecom, Fam110a, Ppia, Atp5c1, Socs2, Rabac1, Paip2, 2900097C17Rik, Coro2a, Ltbp3 
##     Tmem140, Gata2, Pole4, Lamtor5, Gpx4, Bsg, Glul, Ubxn1, Psmb1, Hibadh 
## PC_ 5 
## Positive:  Ypel2, Adgrf5, Dab2, Notch3, Zfp532, Lama3, Pde3a, Wasl, Sema3g, Fkbp10 
##     Kif1c, Jcad, Atp6v0a1, Htra1, Tril, Lrp10, Mecom, Diaph2, Nos3, Piezo1 
##     Kdr, Macf1, Gm5953, Sipa1l2, Plxnd1, Klf3, Errfi1, Sulf1, Hs3st1, Fam8a1 
## Negative:  Oaz1, Clns1a, Cks1b, Tpm3-rs7, Cdc45, Pin1, Sh3bgrl3, Uhrf1, Cops3, Mrps22 
##     Osgep, Unc119, Cenpx, Tpst2, Tubb5, Rps13, Polr2d, Dnph1, Mrps11, Pmf1 
##     Pigu, Incenp, Nme1, Ostf1, Btf3, Rbck1, Lrr1, Nt5dc1, Ube2t, Timm13
DimPlot(object = seu, reduction = "pca", group.by = "orig.ident")

new_order <- c("AEC", "HEC", "T1_pre_HSC", "T2_pre_HSC", "FL_HSC", "Adult_HSC")
seu <- reLevel(x = seu, colable = "orig.ident", level_order = new_order)
DimPlot(object = seu, reduction = "pca", group.by = "orig.ident")

6.2 Cell-type specific TF regulon identification

Here we implemented a function to estimate the TF regulon cell-type specificity based on JS divergence Moran et al. (2010), and then estimated the significance using a bootstrap method.

6.2.1 Get cell-type specific TF regulons for indicated cell type:
cts_tfs_AEC <- findCellTypeSignature(x = atfr, pos_clust = "AEC", clust_name="stage", use_assay ="viper", nperm = 10000, ncores = 2)
cts_tfs_AEC_diff <- cts_tfs_AEC %>% filter(CSS>0.3 & padj <0.05 & exp_percent>0.3)
head(cts_tfs_AEC_diff)

The CSS is the cell-type specific score of each TF represent the cell-type specificity, and the padj determines the significance. Moreover, exp_percent determine the percentage of cells expressed in indicated stage.

6.2.2 Get cell-type specific TF regulon for all cell types:
cts_tfs_all <- findAllCTSRegulons(x = atfr, clust_name="stage", use_assay ="viper", fdr = 1, css=0, nperm = 10000, ncores=2)
cts_tfs_diff <- cts_tfs_all %>% filter(CSS>0.3 & padj <0.05 & exp_percent>0.3)
head(cts_tfs_diff)
6.2.3 Visualization of cell-type specific TF regulons

Here we visualized the cell type specific TF regulons using different methods.

i). Heatmap

library(pheatmap)
act_mat <- assay(x = altExp(x = atfr, e = "viper"))[unique(cts_tfs_diff$regulon),]
anno_col <- as.data.frame(colData(atfr))[,"stage",drop=FALSE]
library(RColorBrewer)
cols <- brewer.pal(3, "YlOrRd")
cols <- c("#1565c0","white","#c62828")
pal <- colorRampPalette(cols)(100)
pheatmap(mat = act_mat,scale = "row",cluster_cols = F, clustering_method = "ward.D2", annotation_col = anno_col, color = pal, fontsize = 5, cutree_rows = 6, show_colnames = F)

ii). RadViz plot

mean_act_score <- summarizeClusterValueRegulon(object = atfr, clust_name = "stage", regulon_name = unique(cts_tfs_diff$regulon), use_assay = "viper")
mean_act_score <- data.frame(mean_act_score,gene_name=row.names(mean_act_score))
plot_frm <- dplyr::left_join(x = cts_tfs_diff[,1:2],y = mean_act_score, by = c("regulon"="gene_name"))
plot_frm$stage <- factor(plot_frm$stage, levels = new_order)
RadvizEnrichPlot(x = plot_frm, anchors = new_order, color_by = "stage", optim_anchor = FALSE, text_size = 4)

6.3 TF regulon-pathway similarity analysis

The TF regulons usually showed functional correlations with pathways. Here we developed a function to performed jaccard test between TF regulon and Reactome pathways.

data(Reactome_pathway)
mm_pathway <- Reactome_pathway$`Mus musculus`
cts_tfs_AEC_regulon <- as.list(Regulon(atfr)[cts_tfs_AEC_diff$regulon])

enrich_res <- calcuModuleSimilar(tfr_list = cts_tfs_AEC_regulon, pathway_list = mm_pathway, ncores = 2)
head(enrich_res)
library(easyalluvial)
library(parcats)
p = alluvial_wide(enrich_res, max_variables = 2, fill_by = "first_variable",
                  colorful_fill_variable_stratum=F,col_vector_value=NA)
parcats(p, marginal_histograms = FALSE, data_input = data, height=500,width = 500)

To simplify the enriched paways for each TF regulon, we clustered the pathway terms based on the similarity of target gene sets, The group column specific the number of groups clustered for each TF regulon, while the firstInGroup indicated the most significant enriched terms in each group.

7. Other parallel analysis

In this section, we give some examples of SCENIC and DoRothEA pipeline.

7.1 Runing SCENIC with simplified steps

We implement a function named runSCENIC() to run SCENIC analysis in a simplified way. First of all, you need to import motif ranking and annotation data from cisTarget databases.

# 1. import data
mm_motifranking <- RcisTarget::importRankings("../ext_data/cisTarget_databases/mm10__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather")
## Using the column 'features' as feature index for the ranking database.
mm_annotation <- data.table::fread("../ext_data/cisTarget_databases/motifs-v8-nr.mgi-m0.001-o0.0.tbl",sep="\t")
# 2. inferring GRNs
atfr <- inferGRNs(x = atfr, method="genie3", use_assay="logcounts", use_regulator = intersect(row.names(atfr),TF_Symbols$mouse), ncores=6)
# 3. runing SCENIC
scenic_auc <- runSCENIC(x = atfr, motif_ranking=mm_motifranking, motif_annotation=mm_annotation, use_assay="logcounts")
## Genes in the gene sets NOT available in the dataset: 
##  Arid5b:     6 (12% of 50)
##  Cebpb:  2 (4% of 50)
##  Creb3:  7 (14% of 50)
##  Deaf1:  1 (2% of 50)
##  Dnmt1:  1 (2% of 50)
##  Elf1:   8 (16% of 50)
##  Gata2:  3 (6% of 50)
##  Gli3:   9 (18% of 50)
##  Hivep1:     3 (6% of 50)
##  Hsf2:   11 (22% of 50)
##  Jund:   6 (12% of 50)
##  Klf12:  9 (18% of 50)
##  Klf3:   4 (8% of 50)
##  Klf5:   7 (14% of 50)
##  Mecom:  6 (12% of 50)
##  Mybl2:  2 (4% of 50)
##  Nfe2l2:     7 (14% of 50)
##  Onecut2:    2 (4% of 50)
##  Patz1:  4 (8% of 50)
##  Pou2f2:     3 (6% of 50)
##  Runx3:  18 (36% of 50)
##  Sall2:  8 (16% of 50)
##  Sall4:  4 (8% of 50)
##  Scmh1:  3 (6% of 50)
##  Stat4:  8 (16% of 50)
##  Stat6:  5 (10% of 50)
##  Tfap4:  8 (16% of 50)
##  Tfdp2:  4 (8% of 50)
##  Tgif2:  9 (18% of 50)
##  Xbp1:   7 (14% of 50)
##  Zbtb22:     13 (26% of 50)
##  Zbtb46:     14 (28% of 50)
##  Zfp105:     4 (8% of 50)
##  Zfp212:     5 (10% of 50)
##  Zfp275:     8 (16% of 50)
##  Zfp408:     14 (28% of 50)
##  Zfp449:     9 (18% of 50)
##  Zfp65:  7 (14% of 50)
##  Zfp780b:    10 (20% of 50)
##  Zfp790:     6 (12% of 50)
##  Zfp831:     11 (22% of 50)
##  Zfpm1:  5 (10% of 50)
##  Zscan20:    9 (18% of 50)
##  Zscan21:    11 (22% of 50)
## Warning in addMotifAnnotation(auc = motifs_AUC, nesThreshold = nesThreshold, :
## The input TFs are not named, all TFs will be used with all Gene Sets.
## Using BiocParallel...
## Quantiles for the number of genes detected by cell: 
## (Non-detected genes are shuffled at the end of the ranking. Keep it in mind when choosing the threshold for calculating the AUC).
##    min     1%     5%    10%    50%   100% 
##  80.00 235.38 286.25 309.30 508.00 750.00
## Warning in .AUCell_calcAUC(geneSets = geneSets, rankings = rankings, nCores =
## nCores, : Ignoring the following empty sets: NA, NA, NA, NA
## Warning in .AUCell_calcAUC(geneSets = geneSets, rankings = rankings, nCores =
## nCores, : Using only the first 50 genes (aucMaxRank) to calculate the AUC.

7.2 Runing DoRothEA

exprMatr <- as.matrix(assay(atfr,"logcounts"))
library(dorothea)
regulons <- dorothea::dorothea_mm %>% filter(tf %in% intersect(row.names(exprMatr),TF_Symbols$mouse))
tf_activities <- run_viper(exprMatr, regulons, options =  list(method = "scale", minsize = 5, eset.filter = FALSE, cores = 1, verbose = TRUE))
## 
## Computing the association scores
## Computing regulons enrichment with aREA